Fluctuation induced vortex pattern and its disordering 
in the fully frustrated XY model on a dice lattice 
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A highly degenerate family of states, in which the plaquettes with the same sign of vorticity 
form clusters of three [proposed in Phys. Rev. B 63, 134503 (2001)] is proven to really minimize 
the Hamiltonian of the fully frustrated XY model on a dice lattice. The harmonic fluctuations 
are demonstrated to be of no consequence for the removal of the accidental degeneracy of these 
states, so a particular vortex pattern can be stabilized only by the anharmonic fluctuations. The 
structure of this pattern is found and the temperature of its disordering due to the proliferation of 
domain walls is estimated. The extreme smallness of the fluctuation induced free energy of domain 
walls leads to the anomalous prominence of the finite-size effects, which prevents the observation 
of vortex-pattern ordering in numerical simulations. In such a situation the loss of phase coherence 
may be related to the dissociation of pairs of fractional vortices with the topological charges ±1/8. 
In a physical situation the magnetic interactions of currents in a Josephson junction array will be 
a more important source for the stabilization of a particular vortex pattern then the anharmonic 
fluctuations. 

PACS numbers: 74.81.Fa, 64.60.Cn, 05.20.-y 



I. INTRODUCTION 

The uniformly frustrated XY model has been intro- 
duced by Teitel and Jayaprakash pj for the descrip- 
tion of a regular array of superconducting islands con- 
nected with each other by Josephson junctions (a Joseph- 
son junction array Q) in the presence of a uniform 
magnetic field. During the last two decades the main 
attention has been concentrated on the investigation 
IlSESSSSHIIJElEIIllIi of so-called fully 
frustrated models (on various lattices), which in terms 
of array correspond to having a half-integer number of 
the superconducting flux quanta per plaquette Q. The 
models belonging to this class can be also used for the 
description of a planar magnet in which the neighboring 
spins can have either ferromagnetic, or antiferromagnetic 
interaction, the number of the antiferromagnetic bonds 
in each plaquette being odd 

The ground states of the uniformly frustrated XY 
models are characterized by the combination of the con- 
tinuous and discrete degeneracies 0. The former is re- 
lated to the invariance of energy with respect to the 
global phase rotation, whereas the latter can be discussed 
in terms of the formation of a particular vortex pattern. 
In the fully frustrated models the numbers of plaquettes 
which contain positive and negative vortices should be 
equal to each other. 

Since vortices of the same signs repel each other, the 
energy is minimized when the vorticities of neighboring 
plaquettes are of the opposite signs. On a square lattice 
this requirement can be simultaneously satisfied for all 
pairs of neighboring plaquettes, which allows one to con- 
clude that the ground state has the checkerboard struc- 
ture and a twofold discrete degeneracy . With increase 
of temperature two different phase transitions can be ex- 
pected to take place , one of which is related to the loss 



of phase coherence and the other can be associated with 
vortex-pattern disordering. The analysis of the mutual in- 
fluence between the two types of topological excitations 
allows one to conclude that in the fully frustrated XY 
model on a square lattice the former has to take place at 
lower temperature than the latter 0] ■ 

Triangular lattice also allows for the construction of a 
doubly degenerate pattern in which the vorticities are of 
the opposite signs for all pairs of neighboring plaquettes 
0, 0| . It turns out possible to show that the thermody- 
namic properties of the fully frustrated XY model on a 
triangular lattice (including the sequence of phase tran- 
sitions) are completely analogous to those of the model 
on a square lattice [R| ■ 

On a honeycomb lattice the situation is more complex, 
because the family of ground states is characterized by 
an infinite accidental degeneracy |6|, which can be de- 
scribed in terms of the formation of parallel zero-energy 
domain walls 9J. In such a case the structure of vortex 
pattern at low, but finite temperatures cannot be deter- 
mined without taking into account the contribution to 
free energy from the small amplitude fluctuations in the 
vicinities of different ground states. This mechanism of 
the removal of an accidental degeneracy 0, is often 
referred to as "order-from disorder". In systems with a 
continuous degeneracy it is usually sufficient to com par e 
the contributions from harmonic fluctuations 

Recently it has been discovered that in the fully frus- 
trated XY model on a honeycomb lattice the order-from- 
disorder mechanism does not work at the harmonic level 
The difference between the free energies of fluctua- 
tions appears only when one takes into account the an- 
harmonicities, and, as a consequence, is proportional not 
to the first, but to the second power of temperature. This 
feature leads to the unusual prominence of the finite-size 
effects E3. 
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The present work is devoted to the investigation of 
the fully frustrated XY model on dice lattice |2l| (see 
Fig. 1). Like square, triangular and honeycomb lattices, 
dice lattice consists of identical plaquettes (which in this 
case are rhombic) and equivalent bonds. Since the in- 
variant description of a Josephson junction array can be 
achieved only in terms of variables which are defined on 
lattice bonds [the gauge-invariant phase differences, see 
Eq. 0}], and not on sites, dice lattice can be consid- 
ered as one of the four basic lattices for the investigation 
of the uniformly frustrated XY models. On more com- 
plex lattices (containing nonequivalent plaquettes) the 
correspondence between an array and a frustrated XY 
model is likely to be broken as a consequence of the phe- 
nomenon of the "hidden incommensurability" , related to 
the redistribution of magnetic field between the plaque- 
ttes by screening currents in asymmetric superconducting 
islands HEI. 

Recently a hypothesis has been but forward ^3 that 
in the ground states of the fully frustrated XY model on 
a dice lattice the vortices of the same sign form three- 
vortex clusters (triads), and a highly degenerate family 
of states has been proposed, which satisfies this crite- 
rion and can be described in terms of the formation of a 
network of intersecting zero-energy domain walls. In Sec. 
II we present a rigorous prove that these states indeed 
correspond to the absolute minimum of energy. 

Sec. Ill is devoted to the analysis of harmonic fluctu- 
ations. We reveal the existence of a hidden gauge sym- 
metry, which allows one to conclude that for a particular 
choice of boundary conditions the set of the eigenvalues 
of the harmonic Hamiltonian is exactly the same for all 
ground states. As a consequence, the free energy of the 
harmonic fluctuations cannot be the source for the se- 
lection of a particular vortex pattern. This conclusion is 
valid also for quantum generalizations of the model. All 
these properties resemble very much the analogous prop- 
erties of the fully frustrated XY model on a honeycomb 
lattice [H|. 

In Sec. IV the lowest order contribution to the free en- 
ergy of anharmonic fluctuations is considered. In partic- 
ular, the gauge symmetry mentioned above is applied to 
demonstrate that the contribution related to the fourth- 
order terms in the Hamiltonian is the same for all ground 




FIG. 1: Dice lattice is the simplest periodic lattice which can 
be constructed from identical rhombic plaquettes with three 
different orientations. 



states and, therefore, is of no consequence for the selec- 
tion of a vortex pattern. The continuous approximation 
is used to show that the fluctuation induced interaction 
of zero-energy domain walls is extremely weak and de- 
cays inversely proportionally to the fifth power of the 
distance between them. The comparison of the numeri- 
cally calculated free energies of anharmonic fluctuations 
in different periodic ground states allows us to establish 
the vortex pattern which can be expected to be stabilized 
at low temperatures, and to find the fluctuation induced 
free energy of zero-energy domain walls. 

In Sec. V we analyze how the vortex pattern selected 
by anharmonic fluctuations becomes disordered when 
one takes into account the fluctuations of another type, 
namely, the formation of domain walls, and propose an 
estimate for the temperature of the phase transition 
which can be associated with the proliferation of such de- 
fects. In this section we also discuss the finite-size effects 
which interfere with the observation of vortex-pattern or- 
dering in finite samples and show that in the considered 
system they are extremely prominent (exactly for the 
same reasons as in the case of a honeycomb lattice). 

In Sec. VI the interplay between the vortex-pattern dis- 
ordering and the loss of phase coherence is considered, 
whereas Sec. VII is devoted to a discussion of another 
mechanism of the removal of an accidental degeneracy 
related to magnetic interactions of currents in the ar- 
ray [23L |24|. In the concluding Sec. VIII our results are 
summarized and compared with the results of numerical 
simulations of Cataudella and Fazio [l4j |. 

The interest to magnetically frustrated systems with 
dice lattice geometry has appeared after Vidal et al. [2{| 
have discovered that at full frustration the ground state 
of a single electron which can jump between the near- 
est sites of a dice lattice is infinitely degenerate and that 
the corresponding wave function can be chosen as an ar- 
bitrary linear combination of an infinite number of ex- 
tremely localized wave functions, each of which covers 
only a finite (and small) number of sites [25|. Since it is 
known that the structure of the superconducting state 
in a wire network is determined (in the mean-field ap- 
proximation) by the structure of the ground-state wave 
function of the single-electron problem with the same 
geometry [26], a suggestion has been put forward that 
at full frustrateion the superconducting state in a dice 
network may have a disordered (glassy-like) structure 
|27| . However, recently it has been shown [24| that 
the inclusion into analysis of the forth-order term of the 
Ginzburg-Landau functional strongly decreases the am- 
biguity in the determination of the structure of super- 
conducting state in a fully frustrated wire network with 
the dice lattice geometry. The set of states minimizing 
the free energy of such a network turns out to be in one- 
to-one correspondence with the set of the ground states 
of the fully frustrated XY model discussed on this arti- 
cle. In recent years magnetically frustrated wire networks 
and Josephson junction arrays with the dice lattice ge- 
ometry have both been the subject of active experimental 



3 



investigations pi li^ . 

II. THE MODEL AND THE GROUND STATES 

A. The definition of the model 

A uniformly frustrated XY model [lj can be defined 
by the Hamiltonian [3(j 

H = - cos(ip k - (fj - A jk ) , (I) 
(jk) 

where the summation is performed over all bonds (jk) of 
a regular two-dimensional lattice. In terms of a Josephson 
junction array J is the Josephson coupling constant of a 
single junction, fluctuating variables fj are the phases of 
the order parameter on superconducting grains j forming 
the array, whereas quenched variables, 

A jk = — rfrA(r), (2) 

arc defined by the integral of the vector potential A(r) of 
the external magnetic field along the bond (jk), <3>o being 
the superconducting flux quantum. The form of Eq. 
assumes that the currents in the array are sufficiently 
small, so their proper magnetic fields can be neglected. 

When the magnitude of the field corresponds to a half- 
integer number of flux quanta per plaquette, the directed 
sum of Ajk = —Akj along the perimeter of a plaquette 
in positive direction (which below is designated as J2n) 
has to satisfy the constraint 

A ]k = 7T (mod 2tt) (3) 

□ 

on each plaquette of the lattice. In such a case the model 
is called fully frustrated 3]. In a more general case of 
a uniformly frustrated XY model, the right hand-side 
of Eq. (PJ) should be replaced by 27r/(mod 27r), where 
the frustration parameter / describes the magnitude of 
the external magnetic field in terms of the number of 
flux quanta per plaquette. It is sufficient to consider the 
interval / € [0, i], because all other values of / can be 
reduced to this interval by a simple replacement of vari- 
ables 0. The term "fully frustrated" is used for the case 
of / = 1/2, the maximal irreducible value of /. 

Since both variables ipj and variables Ajk depend on 
a choice of a gauge, it is often more convenient to de- 
scribe different states of the system in terms of the gauge- 
invariant phase differences, 

Ojk = fk - fj - Ajk = ~0 k j , (4) 

defined on lattice bonds. Below we will always assume 8jk 
to be reduced to the interval (— w, n). It follows from the 
definition of these variables that in the fully frustrated 
model they have to satisfy the constraints, 

£ 9 3 k = 7T (mod 2tt) , (5) 
□ 



completely analogous to Eq. J3J). One usually says that 
a given plaquette contains a positive (or negative) half- 
vortex when the left-hand side of Eq. © is equal to +n 
(or — 7r). Different minima of \Q (including the ground 
states) can be then identified in terms of a corresponding 
vortex configuration. 

The variation of Eq. with respect to fj results in 
the current conservation equation for the site j, the value 
of the current in the junction (jk) being given by 

Ijk = Iasmdjk , (6) 

where Iq = (2e/h)J is the critical current of a single 
junction. 



B. Minimization of energy 

Dice lattice (Fig. I) is formed by two types of sites, 
with the coordination numbers three and six, connected 
with each other in such a way that each bond connects 
two sites with different coordination numbers. Below we 
will always use index k to denote the threefold coordi- 
nated sites of a dice lattice and index j to denote the six- 
fold coordinated sites. For example, the bond (jk) con- 
nects the sixfold coordinated site j with the threefold 
coordinated site k. 

The minimization of the Hamiltonian with respect 
to all variables ip k for the given values of the variables fj 
can be performed exactly. To describe the result of this 
procedure it is convenient to introduce also the gauge- 
invariant phase differences Xjj' = ~Xj'j defined on the 
bonds of the triangular lattice T formed by the sixfold 
coordinated sites j. A natural way to do it consists in 
requiring that for each triangle formed by the sites j, j' 
and k (where j and j' are the nearest neighbors of k) the 
sum of the three gauge-invariant phase differences taken 
along its perimeter in the positive direction should be 
equal to ±7r/2 modulo 2ir, where the sign should be the 
same for all triangles. In what follows we assume this sign 
to be negative, 

Xjj' + Qj'k + Okj = -tt/2 (mod 2tt) . (7) 

Since Xjj' — ~Xj'ji t ne constraint (jSJ) then automatically 
follows from Eq. Q. On the other hand, on each triangu- 
lar plaquette of T the directed sum of the variables Xjj' 
has to satisfy the constraint 

$3^,= tt/2 (mod 27r), (8) 

□ 

which can be obtained by summation of Eq. Q for three 
neighboring triangles with the common cite k. 
The minimization of 

3 

Ek = -J^2cos9j ak , 

a=l 
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(where j a with a — 1, 2, 3 are the three nearest neighbors 
of k on a dice lattice numbered in the positive direction) 
with respect to (fk for the given values of Xjih > Xjaja ancl 
Xi 3 ji satisfying the constraint iJSJ gives 



where 



Y, 



E k = -Jy/3 - 2Y k 



l Xjij2 + siD -Xj 2 j 3 + sin Xj 3 ji < 3 / 2 



Since E(Y) = —Jy/3 — 2Y is a concave function of Y 
and the sum of the variables Y k over the whole lattice 
with the periodic boundary conditions should be equal 
to zero, the absolute minimum of 



H=Y J E{Yk) 



on such a lattice is achieved when Y k = for all k. In 
the next subsection we demonstrate that this require- 
ment can be simultaneously satisfied on all plaquettes. 
Accordingly, the value of energy in the absolute mini- 
mum is given by E = -2y/3JN, where N is the total 
number of the sixfold coordinated sites. 



C. Construction of ground states 

In the case of an isolated triangle the system of two 
equations 

X1+X2+X3 = tt/2 , 
sinxi + sinx2 + sinx3 = 0, 

for three variables Xl> X2 and X3 nas an infinite num- 
ber of solutions. However, the requirement to match the 
solutions on all triangular plaquettes of T leads to the 
removal of this continuous degeneracy. 

Since the form of the constraint JSJ corresponds to hav- 
ing one quater of flux quantum per plaquette, the min- 
imal elementary cell whose periodic repetition allows to 
construct a periodic solution consists of 4 triangles. Such 
a solution can be described by 6 variables \i [defined, for 
example, like is shown in Fig. 2(a)], which have to satisfy 
3 constraints of the form (JSJ) : 



Xi +X2 
X4 + X5 
"X3 - X4 



X3 
X6 
X5 



tt/2, 
tt/2, 
tt/2, 



and 3 equations of the form Y k = 0: 



smxi 
sin X4 
sinx3 



smx2 
sin X5 
sin X4 



The fourth constraint, 

-Xi - X2 



X6 



sinx3 = 
sinxe = 
sinxs = 

-3tt/2 , 



(9a) 
(9b) 
(9c) 



(10a) 
(10b) 
(10c) 



and the fourth equation, 

- sinxi - sinx2 - sinxe = . 

are then satisfied automatically, as follows from the sum- 
mation of Eqs. 10 and Eqs. (|10|) . respectively. 

We have found that the numerical solution of the sys- 
tem of six equations by an iterative method al- 
ways converges to the solution shown in Fig. 2(b) or to 
other analogous solutions in which variables x are equal 
to 0, — 7r/4 and 37r/4 on one half of triangular plaque- 
ttes and to 7r, 7r/4 and 37r/4 on the other half. Note that 
each variable x belongs to two neigboring plaquettes, but 
manifests itself on them with the opposite signs. The 
same results are also obtained when one assumes that 
the elementary cell has a different shape shown in Fig. 
2(c). 

Quite remarkably, an attempt to construct a solution 
with a larger elementary cell leads to an overdefined sys- 
tem of equations. For example, an elementary cell con- 
sisting of 8 triangles requires to introduce 12 variables 
Xi which have to satisfy 7 independent constraints of 
the form JSJ) and 7 independent equations of the form 
Yfe = 0. Apparently, one cannot expect a system of 14 
equations for 12 variables to have a nontrivial solution 
which cannot be constructed from the solutions obtained 
for a 4-triangle elementary cell. Thus, the reduction of the 
problem to a triangular lattice has allowed us to make a 
conclusion on the size of the elementary cell which would 
be hardly possible in the framework of analysis in terms 
of the original variables 9jk defined on the bonds of a dice 
lattice. 

In Fig. 2(d) the structure of the elementary cell of Fig. 
2(b) is shown in terms of the variables 9j k . Here sin- 
gle, double and triple arrows correspond, respectively, to 



X3 X2 X6 X5 X3 



(c) 



A 

X3 X5 



3 »/r /4 7r /4 7 /-»\ 

3^/4 7T 3tt/4 A A 

/ \/ V X3 X2 X6 X5 

W4^W4^ J \J \ 

Xi » v < Xi— ^ 




FIG. 2: Construction of a periodic ground state: (a) a possible 
structure of an elementary cell; (b) one of the solutions of 
Eqs. (c) an alternative elementary cell; (d) the same 

elementary cell as in (b), but in terms of 9jk- 
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three different values of djk, 

I 1 1 

7i 3 = arccos | — = ± — = 



1 



satisfying the constraints 

2 -0 1 = 7r /4, 0i+03=7r/2, 2 + 03 = 3tt/4 , 

which lead to the fulfilment of Eq. © on all rhombic 
plaquettes, and the current conservation equation, 

sin 0i + sin 02 = sin 03 , 

the form of which follows from Eq. 10. 

The most compact way of illustrating a structure of 
a given state (a local or global minimum of the Hamil- 
tonian) consists in showing which plaquettes are occu- 
pied by positive and which plaquettes by negative half- 
vortices. In Fig. 3(a) this approach is used to demonstrate 
the structure of the ground state which is obtained by the 
periodic repetition of the elementary cell shown in Fig. 
2(d). Notice that positive and negative half- vortices are 
grouped into clusters of three (triads). The rules which 
allow one to restore the values of 9jk for each bond from 
the structure of vortex pattern (in a ground state) can 
be found in Ref . ll2L 

Since positive and negative half- vortices can be consid- 
ered as occuping the cites of the dual lattice (which in the 
present case is a kagome lattice) , the structure of a given 
state of the fully frustrated XY model on a dice lattice 




FIG. 3: The structure of some ground states. The plaquettes 
with positive vorticities are marked in black. 



can be compared with the structures of different states 
of the antiferromagnetic Ising model on a kagome lattice. 
In particular, according to Wolf and Schotte [sHi hi the 
framework of the Ising model the state with the structure 
shown in Fig. 3(a) is selected when Ji 3> J 2 ^S> J 4 > 
and all other couplings are equal to zero. Here Ji is the 
coupling constant for ith neighbors on a kagome lattice, 
and we have used the notation of Ref. |2J for the clas- 
sification of dice lattice plaquettes as neighbors of each 
other. 

The ground state whose structure is shown in Figs. 
2(d) and 3(a) allows for creation of infinite domain walls 
which brake the periodicity of this state but does not cost 
any energy |l2j| . These zero-energy domain walls can be 
of the two types. 

Fig. 3(b) shows an example of zero-energy domain wall 
of the type I. Note that the orientations of triads formed 
by negative half-vortices (white plaquettes) are different 
above and below the wall. The configuration of arrows 
(defining the values of the variables 9jk) after crossing 
such a wall can be obtained from the old configuration 
(in the absence of the wall) by its reflection with respect 
to a line which is perpendicular to the wall and subse- 
quent inversion of all arrows. There can exist an arbitrary 
number of such domain walls in parallel to each other. By 
creating them at every possible position one obtains an- 
other periodic ground state shown in Fig. 3(c). 

Fig. 3(d) shows an example of zero-energy domain wall 
of the type II. After crossing such a wall the variables 
are changed (in comparison with what they would be in 
the absence of the wall) according to even more simple 
rule, which can be codified as 

01 - * 03 7 03 - * 01 7 S-2 — * —02 ■ (11) 

Note that both black and white triads have different ori- 
entations on two sides of the wall, whereas at the wall 
the shape of white triads is modified. 

Again, there can exist an arbitrary number of the type 
II domain walls parallel to each other. By inserting them 
at every possible position one obtains one more periodic 
ground state shown in Fig. 3(e), in which all triads have 
the modified shape. Alternatively, one can start the whole 
construction from the periodic state of Fig. 3(e) and ob- 
tain the state of Fig. 3(a) by introducing a dense sequence 
of domain walls on crossing which the same rule, Eq. Ullfl . 
is applicable. 

The zero-energy domain walls of different types can 
cross each other without increasing energy. However, it 
follows from the rule for the transformation of the state 
induced by the type I domain wall (described above) that 
a type II domain wall should change its orientation by 
7r/3 each time it crosses a type I domain wall [see Fig. 
3(f)]. A dense network of zero-energy domain walls of 
both types constructed on the background of state (a) 
leads to the periodic ground state shown in Fig. 3(g). 
The structures of the periodic states (a), (c), (e) and (g) 
in terms of the variables Ojk are shown in Fig. 3 of Ref. 
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Note that the formation of zero-energy domain walls is 
related to the changes of the orientation of vortex triads 
(and, in the case of type II domain walls, also of their 
shapes), but does not lead to the appearance of vortex 
clusters of other sizes. 



Eq. i|14b[) into Eq. (|14a() and can be written as 



Any- 



where 



E 

j'=j'(j) 



(15) 



III. HARMONIC APPROXIMATION 

A. Two families of eigenmodes 

The Hamiltonian describing the harmonic fluctuations 
in the vicinity of one of the ground states described in 
Sec. II. C can be written as 



A(w) = (2M 3 + M 6 )u 2 - M^± u * , (16) 



Js 

j' (j) are the six nearest neighbors of j on T and 

Kjj' = (Jjk'Jfk' + Jjk"Jj'k")/Js , 



(17) 



(12) 



Uk) 



k' and k" being the two threefold coordinated sites be- 
longing to the same rhombic plaquette as j and j'. 

The right-hand side of Eq. (|15f) has exactly the same 
form as in the equation describing harmonic fluctuations 
on a triangular lattice with the nearest neighbor interac- 



where Uj are deviations of the variables & from their tion characterized by the coupling constants K jr defined 



equilibrium values on sixfold coordinated sites, v k are 
analogous deviations on threefold coordinated sites and 
the coupling constants Jjk = J cos 9j k acquire one of the 
three possible values 



J , J 2 = ^ 



J_ 

71' 



(13) 



in accordance with the value of 9jk on the bond (jk). 

If phase dynamics in a Josephson junction array can 
be assumed to be non-dissipative and the capacitance 
matrix of the array has only diagonal elements, the lin- 
earized equations of motion for the variables Vk and Uj 
following from Eq. (|12ll can be written as 

(2J S - Mquj 2 )^ = J JkVk , (14a) 

k=k(j) 

{J s - M 3 u 2 )v k = J Jk u j , (!4b) 

j=j(k) 



by Eq. l(T7|l . Quite remarkably, in all the ground states de- 
scribed above these coupling constants acquire only two 
values 



K 2 



2JiJ 3 



J 



Js 3V3 ' 
(Jj + J 3 ) J 2 2J 



Js 



3V3 



(18a) 
(18b) 



which differ from each other by the factor of 2. When a 
half-vortex in the plaquette (jk'j'k") is the central vor- 
tex of a triad to which it belongs, Kjj> = K\, whereas 
otherwise Kjy = K 2 - In all the ground states described 
in Sec. II. C these couplings are distributed between the 
bonds of T in such a way that in each triangular pla- 
quettes one bond has Kjji = K%, whereas the two other 



bonds have K 



jj' 



Kn. 



B. Comparison of different ground states 



where Js = Ji + J2 + J3, M t = (h/2e) 2 d (where i = 
3, 6), and j(k) denotes the nearest neighbors of k. In Eqs. 
(|14|l we have performed the Fourrier transformation to 
the frequency representation and have assumed that the 
self-capacitances of the superconducting islands, C3 and 
Cq, are different for the two types of islands. 

Note that in all ground states discussed in Sec. II. C the 
coupling constants Jjk always have the same three values 
( Ji, J2 and J3) on the three bonds (j a k) connected to any 
site k, as a consequence of which the coefficient standing 
in the left-hand side of Eq. (|14b() does not depend on k. 
This allows to conclude that all eigenmodes with Uj = 
should have the same eigenfrequency lu — (Js/M^) 1 / 2 . 
In the thermodynamic limit the degeneracy of this eigen- 
frequency is equal to one third of the total number of 
modes. 

The spectrum of the brunch with Uj 7^ can be found 
from the equation which is obtained after substitution of 



In the state (e) all sixfold coordinated sites have the 
identical environment in terms of the coupling constants 
Jjk.. As a consequence, the values of the coupling con- 
stants Kjji in this state depend only on the orientation 
of the bonds (jj')- For two of the three possible orienta- 
tions of the bonds they are equal to each other, see Fig. 
4(a). If all sixfold coordinated sites j are renumbered by 
pairs of integers with the same parity (n, m) as shown in 
Fig. 5, Eq. 1151) in this state can be rewritten as 



[2K S - A(cj)]u„, m = 

,(m) 



(m) 



K. 



= K{ 'Un+i, m+ i + K. 
js-{m—l) 

+ A' 2 [u n+ 2,m + Un-2,m] , 

where Ks = K\ + 2K 2 and 



n— l,m-t-l 



(m-1) 



"n+l.m— 1 



(19) 



K 



(to) 



K 1} 4 m) = K 2 



(20) 
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for all values of m. 

For given by Eq. Eq. (JTSJl can be easily 

solved after performing the Fourrier transformation. Sub- 
stitution of u n ^ m oc expi(qn + pm) into Eq. I|19f) gives the 
dispersion relation in the form 



A(lu) = 2K S - 2Ki cos(g + p) 

— 2K 2 [cos(g — p) + cos 2q] 



(21) 



Note that Eq. I|21|l is of the second order in u 2 , which 
corresponds to the existence of the two momentum- 
dependent eigenfrequencies, u>x(q,p) and u>2(q,p), for 
each point in the Brillouin zone in addition to 
w o = (Js/Ma) 1 ' 2 discussed above. 

According to Eq. iJHJ, each zero-energy domain wall 
of the type II leads to the permutation of the coupling- 
constants J\ and J3 in the Hamiltonian of harmonic fluc- 
tuations. Such a permutation does not change the factor 
Js — M3UJ 2 in the right-hand side of Eq. 114b(l , and there- 
fore does not change neither the degeneracy nor the fre- 
quency, ujq, of the family of the eigenmodes with u = 0. 
Since Eqs. Ijl8|l are also invariant with respect to the per- 
mutation of J\ and J3, such a permutation introduces 
no changes to the form of Eq. (|15fl as well. This means 
that the whole set of the eigenfrequencies of the harmonic 
fluctuations in the system does not feel the presence of 
the type II domain walls and is exactly the same in all 
the states which can be obtained from each other by the 
insertion of the type II domain walls. For example, the 
dispersion relation (|21|l is valid not only in the state (e), 
but also in the state (a) , which is characterized by exactly 
the same pattern of the coupling constants Kjji shown 
in Fig. 4(a). 

The zero-energy domain walls of the type I lead to a 
more complex permutation of coupling constants. How- 
ever, in the terms of the coupling constants Kjji the con- 
sequences of this permutation are rather simple and re- 
duce to the permutation of K\ and K 2 for those two 




FIG. 4: The distribution of coupling constants Kjji between 
the bonds of T in different states: (a) states a and e; (b) 
a single type I domain wall; (c) states c and g. Thin lines 
correspond to Kjji = K\ and thick lines to Kjji = -Re- 



orientations of the bonds (7V) that are not parallel to 
the direction of the wall [^3], see Fig. 4(b). This means 
that Eqs. (120ft should be valid only when d m , the number 
of the type I domain walls situated at w! < m, is even 
and should be replaced by 



K 



(m) 



Ko , K. 



(m) 



when d m is odd. 

In the presence of periodic boundary conditions in the 
horizontal direction and open boundary conditions in the 
perpendicular (vertical) direction the irrelevance of such 
permutations of coupling constants for the set of eigen- 
frequencies can be easily demonstrated in the framework 
of the mixed representation which is obtained after per- 
forming the Fourrier transformation with respect to the 
variable n, keeping the variable m as it is. In the terms 
of the variable u m (q) defined in such a way, Eq. Ijl9(l can 
be rewritten as 

[2K S ~ A(co)]u m (q) = K^ m \q)u m+1 {q) (22) 
+ 2K 2 cos(2q)u m (q) + K^-V^Um^iq) , 

where 

Kt m \q) = K (q)exp[i<*{q)8 m ] . 

The dependence of K {m \q) on m enters only through 
s m = (— l) dm = ±1, whereas both 



and 



K (q) = \Ki exp(iq) + K 2 exp(-i<?)| 



a(q) = arg[i^i exp(ig) + K 2 exp(— iq)] 



are independent of m. 

With the help of the simple gauge transformation 



u m {q) = cxp 



«(«) E 



u' m {q) , (23) 



which cannot not change the eigenvalues, Eq. (|22|) can 
be transformed to the form 



[2K S - A(uj)]u' m {q) = K (q)u' m+1 (q) 

+ 2K 2 cos(2q)u' m (q) + Xb(?)«m-i(?) - 



111 

4 



(24) 



2 




1 



6 



FIG. 5: The numbering of the sixfold coordinated sites by 
pairs of integers (n, m) with the same parity. 
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in which all coefficients do not depend on m. This prop- 
erty of Eq. (|2*4"|) proves that for the considered boundary 
conditions the whole set of eigenfrequencies is insensitive 
to the presence of domain walls of the type I. Since Eq. 
(|24|l has been derived from Eq. I|19f). the form of which 
does not depend on the presence of the zero-energy do- 
main walls of the type II, the same conclusion is appli- 
cable also in the presence of an intersecting network of 
zero-energy domain walls of both types. 

Note that this does not mean that the spectrum of 
fluctuations in an infinite system, understood as the de- 
pendence of the eigenfrequencies on the two-dimensional 
wave- vector, is the same for all periodic ground states. On 
the contrary, it should have different forms in the states 
whose transformation into each other requires the inser- 
tion of a periodic sequence of the type I domain walls. 
In particular, the dispersion relation in the state (c) and 
state (g) [which can be obtained from each other by the 
insertion of the dense sequence of type II domain walls 
and both are characterized by the distribution of the cou- 
pling constants Kjj, shown in Fig. 4(c)] can be obtained 
by substitution of u m (q) oc exp(ipm) into Eq. (|24|l and 
is of the form 



A(w) = 2K S - 2K (q) cosp - 2K 2 cos 2q , 



(25) 



Apparently, it does not coincide with the dispersion re- 
lation in the states (a) and (e) given by Eq. (|21() . Note 
that the elementary cell of the state (c) consists of twelve 
sites, and, therefore, a more straightforward approach to 
the derivation of its dispersion relation would give it in 
the form of the determinant of a twelve by twelve matrix. 

However, the free energy of harmonic fluctuations, F2, 
which in a general situation (that is, when the quantum 
effects are also taken into account) can be written as 



r£ln(2 S inh^ 



(26) 



is determined entirely by the set of the eigenfrequencies 
of the system {to}. Thus, our results demonstrate that 
for the considered boundary conditions the value of F2 is 
exactly the same for all the ground states discussed above 
even in the case of a finite system. For other types of 
boundary conditions the same property will be recovered 
in the thermodynamic limit. 

Naturally, these conclusions should remain valid both 
in the zero-temperature limit, when Eq. (|26|l is trans- 
formed into the expression for the energy of zero-point 
fluctuations, 



E 2 = - 



dq dp 
(2^P 



UJQ + LJl (q,p) + Ul>2 (<?, p)] 



and in the classical limit (h — > 0), when the dynamical 
properties of the system are of no importance, and Eq. 
(|26l) is reduced to 



Fo = - 



1/ 


f dqdp 






J W ^ 


T 



where A(q,p) denotes the function of q and p standing 
in the right-hand side of the corresponding dispersion 
relation, Eq. J2H) or Ea.pB). 

Another system, in which the accidental degeneracy 
of its ground states remains unbroken when the free en- 
ergy of the harmonic fluctuations is taken into account, is 
the fully frustrated XY model with a honeycomb lattice 
[15j. The method used in this section (the construction 
of the gauge transformation which reduces the linearized 
equations of motion for fluctuations in different ground 
states to the same form) presents a generalization of the 
approach of Ref. 0, where the analogous gauge trans- 
formation has been constructed for the harmonic part of 
the Hamiltonian. 

The analysis of this section can be generalized for the 
case when the capacitance matrix of a Josephson junc- 
tion array in addition to the self-capacitances of super- 
conducting islands also takes into account C\ , the capac- 
itances of Josephson junctions forming the array. This 
will lead to the replacement 

J a -> J a (u) = J a ~ M iu 2 , 

where Mi = (fr/2e) 2 Ci, in all dispersion relations, but 
will not bring any changes in the qualitative conclu- 
sions. It is equally possible to include into considera- 
tion an Ohmic dissipation of each junction described by 
a frequency-dependent harmonic contribution to its Eu- 
clidean action. 



Correlation functions 



Since the evolution of the variables Uj can be described 
by Eqs. I|15|l. which contain only coupling constants Kjj, , 
but not the original coupling constants Jjk, the symme- 
try of the correlation function of the variables Uj in a 
given state will be determined entirely by the symmetry 
of the configuration of Kjj, in this state. In particular, 
in situation when the value of Kjj, depends only on the 
orientation of the bond (Jj') and can be equal only to 
K\ or K2, the correlation function 



G, 



((Uj - u r ) 2 ) 



(27) 



where j and j' are the nearest neighbors of each other on 
T, should also be dependent only on the orientation of 
the bond (jj') and acquire one of the two possible values, 
which below are denoted G\ and G2. The same property 
can described by saying that in such states Gjj> depends 
only on whether Kjj, is equal to K\ or to K2, 



G 



jj' 



Gi for Kjj, = K x , 
G 2 for K jjt = K 2 . 



(28) 



The type II domain walls are simply of no consequence 
for the values of Kjj, and, therefore, for the correlation 
functions of the variables u. On the other hand, it is 
rather evident that the gauge transformation (|23|l leaves 
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invariant the form of the correlation functions of the vari- 
ables u n>m with the same value of m. This means that 
the correlation function ((uj — Uji) 2 ) cannot be sensitive 
to the presence of the type I domain walls which do not 
pass between the points j and f . Since the sites j and f, 
which are the nearest neighbors of each other on T, are 
too close to have a domain wall passing between them, 
Eq. i|28[l will be applicable also in the presence of an ar- 
bitrary set of zero-energy domain walls. 



IV. ANHARMONIC FLUCTUATIONS 

A. Basic relations 

In this section and below our analysis is restricted to 
the classical version of the model. It is well known that in 
a classical system the leading contribution to the free en- 
ergy related to anharmonic fluctuations, -F an h = F3 + F4,, 
is the sum of the two terms, 



F 3 = —. 



1 

2T 



H (3) 



and 



F 4 = (ff (4) ) 



(29) 



(30) 



where and are, respectively, the third- and the 
forth-order contributions to the expansion of the Hamil- 
tonian in the vicinity of a particular ground state, and an- 
gular brackets denote the averages over thermodynamic 
fluctuations calculated with the help of the harmonic 
Hamiltonian. 

In the considered problem H ^ and can be writ- 
ten as 



ff (3) = W(3). 



k 



= -i 



J 'ik( U 3 ~ V kY 



and 



(4) 



24 



Jjk(uj-v k ) 4 



where coupling constants Jjk are exactly the same as in 
the harmonic Hamiltonian, Eq. whereas coupling 

'jk acquire one of the six possible 
±Jsin# a in accordance with the value of 
vjk on the bond (jk). Due to the current conservation 
condition coupling constants Jj fc have to satisfy the con- 
straints, 



constants Jj fe = Jsin 
values Jj k = 



£ J' jk - , £ J' jk - , 

j=j(k) k=k(j) 



(31) 



on all sites of the lattice. 



B. Invariance of F4 

In the case of a dice lattice there exists a convenient 
way to separate the fluctuations on the two types of sites 
from each other, which allows to considerably simplify 
the calculation of the averages in Eqs. ffify and (|3Tjjl . It 
consists in the replacement of variables 



, (0) 
Vk = w k + v k 



= J 3k U j/J S , 

3=j(k) 



(32) 



which transforms the harmonic Hamiltonian (|12H into 



\ JsW k 



(33) 



Uj') 



where the coupling constants Kjji are, naturally, the 
same as have been obtained in Sec. III.B, see Eqs. <|17[) 
and (|18|l . after the exclusion of the variables v k from the 
equations of motion. 

A simple form of Eq. (|33|) , in which each variable Wk is 
decoupled from all other variables, makes the calculation 
of averages with respect to the fluctuations of Wk a very 
straightforward procedure. Substitution of Eq. (|3^|) into 

the expression for with subsequent expansion of the 
result in powers of Wk allows one to express the result of 
such an averaging of as a fourth-order polynomial 
of the variables 



Pi { u jn u j 2 ) u j3 ) 



32 J ""33 I 

-1^J Q (3< W 2 ) 2 + 6< W 2 )^+^ 

a=l 



(34) 



where (w 2 ) — T / Jg is the value of 
same for all k, whereas 



. w k) 



which is the 



,(°) 



■Is 



All terms which are odd in Wk have disappeared from Eq. 
(|34|l due to the corresponding symmetry of the Hamilto- 
nian l|33|l . Note that the form of ^4(^1, ^2,1*3) does not 
depend on fc. To achieve that we have renumbered in Eq. 
(|34|l the three sites j which are the nearest neighbors of 
k as j a = j a (k) in such a way that Jj a k = Ja- 

Since Hamiltonian (|33|l does not include linear terms, 
the result of the Gaussian averaging of P4 (uj 1 , Uj 2 , Uj 3 ) 
will have a form of a second order polynomial, 
P2(Gj 1 j 2 ,Gj 2 j a ,Gj 1 j 3 ), whose three arguments are the 
nearest neighbor correlation functions defined by Eq. 
(|27|l . Instead of looking for the explicit form of P2, it 
is sufficient to notice that since the central vortex of 
a triad is always surrounded by the bonds with Jjk 
equal to Ji or J3, we will always have Kj ± j 3 — K\ and 
K3131 = K3233 — and, as a consequence of Eq. JSSJl, 
the result of the averaging of will have the same 
form, 

r(4) 
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for all k independently of what particular ground state is 

considered. Therefore, the value of F 4 = J2k(H^) will 
be exactly the same for all the ground states which we 
are trying to compare already at the level of the separate 
terms in this sum. 



C. Simplification of F3 

(3) 

The result of the averaging of HI with respect to 
fluctuations of Wk can be in analogous way reduced to 
the sum of two terms, the first of which, 



P 1 (u jl ,u j2 



U 3:J = 



3 

' E J lk U 0o 
a=l 



is a first-order and the second, 



1 3 

gE J ']ak u 3 
a—1 



a third-order polynomial of the variables Uj a , which are 
assumed here to be numbered in the same way as in Eq. 
I)34p. Both Pi and P 3 depend on k only through the fac- 
tor Tfc = ±1, which, for example, can be chosen to be 
determined by the sign of Jj lfc . 

The sum of Pi (uj 1 , Uj 2 , uj 3 ) over k can be reordered as 
a sum over j, 

„2\ 



fc 



Pi( Ujl ,u 3 



{w 



where 



k=k(j) 



[J'jk(Jj'k + Jj"k) - Jjk(Jj> k + Jj"k)} 



and j' and j" are the other two nearest neighbors of k in 
addition to j. It is not hard to notice that all coefficients 
Cj are equal to zero as a consequence of Eqs. J3TJ. 
This allows one to express the result of the averaging 

of [i?( 3 )] over fluctuations of the variables Wk as 



H 



(3) 



= R 2 



E^(< 



'31 1 u 32 I 



where 



R = ^2P 3 (u jl ,u h ,u h ) 



(35) 



(36) 



and Pq is a sixth-order polynomial of its arguments, the 
form of which is the same for all k. Exactly like it happens 
with P4, the averaging of Pq (uj 1 ,Uj 2 ,Uj 3 ) over fluctua- 
tions of Uj produces the same expression for all k inde- 
pendently of what particular ground state is considered. 
That means that any difference between the values of 
F an h in different ground states can result only from the 
first term in Eq. (|35|l . 



-Fanh = const 



— (R 2 ) 



(37) 



D. Explicit expressions 

Let us start with comparing the free energies of anhar- 
monic fluctuations in the states (a) and (e). Instead of 
calculating F an h separately for both these states it is more 
convenient to construct an explicit expression directly for 
^-^Tnh = ^anh _ ^anh' tne difference in i^nh between the 
states (e) and (a). Since both these states are character- 
ized by the same form of the effective Hamiltonian for the 
variables uj , the construction of such an expression does 
not require the application of the gauge transformation 
given by Eq. I|23l) . The same is true also for the whole 
set of states which can be constructed from the state (a) 
or state (e) by the insertion of some sequence of type II 
domain walls. 

If the state (e) shown in Fig. 3(e) is rotated by 7r/3 in 
such a way that the orientation of the possible type II 
domain walls becomes horizontal, the expression for R, 
Eq. I|3t)|) . in this state can be rewritten as 



Re 



(38) 



where we have again used the numbering of sites defined 
by Fig. 5, a m = (-ir, 



n=m (mod 2) 

"f" { U n-\-l,7n-\-l-> ^71,171-, U n — 1,771+0] ) 

and = TfcP3 is an invariant version of P3, 



(39) 



D 3 + (U1,U2,U 3 ) 



1 3 

6Jl^ J ' a 

a=l 



E J b( u a ~ U b ) 



(40) 

the form of which does not depend on k. Instead of intro- 
ducing a definition simply for S^, we have used Eq. Ij39(l 
to define a more general object S/^, where superscript /i 
can be equal to ±1 or 0. However, for the compactness of 
notation we will usually replace \x — ±1 simply by plus 
or minus. In Eq. Ij40|l the constants J' jk are expressed in 
terms of three constants, 



J! 2 = J sin 2 = -|J, 



Jo = — J sin ( 



1 1 

V3 V6 



(41a) 
(41b) 
J, (41c) 



whose sum is equal to zero in accordance with Eqs. I|31|l . 

Substitution of Eqs. lO and Eqs. (JUI into Eq. l(4T)|) 
allows to reduce it to 

P^(ui,u 2 ,u 3 ) = K 3 (u 2 - u 3 ) 2 (u 3 - ui) (42) 
+ K' 3 (ui - u 2 f + K' 3 \u 3 - mf , 
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where 



A : 



6%/rj' Ka 9%/6 



3 K" ( 1 



V9V3 36 V6 



However, in all the ground states that we consider the last 
term from Eq. (|42|l after substitution of P 3 = TkP 3 into 
Eq. I|36|) cancels with analogous term from the neighbor- 
ing triangular plaquette, which allows one to put K 3 = 0. 

Each time a type II domain wall is crossed the replace- 
ment of variables 6j described by Eqs. Hll|) has to take 
place. In terms of the expression for R this procedure is 
translated into the replacement of each term of the form 
P 3 + («i, u 2 ,u 3 ) by 

Pf(u!,U 2 ,U 3 ) = --P 3 + ("3,M2,Ml) ■ 

As a consequence, the value of R for the state which is 
obtained after the insertion of the dense sequence of type 
II domain walls (the structure of this state is obtained 
after rotating Fig. 3(a) by 7r/3) will have the form 



(43) 



Substitution of Eqs. (JSHJ and ((EJ into Eq. allows 
to express SF^ as 



E. Continuous approximation 

When Eq. (|46|) is substituted in the expression for S 1 ^ 
given by Eq. (|3*9"jl with fx — 0, all terms proportional to 
K3 cancel each other after the summation over n, so it 
becomes possible to rewrite this expression as 



'->,,, 



4i<T, 



E (v««) 2 (v; 



I m' =m+l 



n—ra' (mod 2) 



E ( v « u ) 2 ( v 

n—m' (mod 2) 



m / | rn'—m 



(49) 



where 



is a lattice analog of the derivative in the horizontal di- 
rection and 



± 



U 



n+l,m'=pl 



^aTh^^E^- 1 ) 



1=1 



where V(m)) is the average 

2 



V(m) 



\ tJ rni ni2 / ' 



(44) 



(45) 



which depends only on to = toi — TO2, is the size of 
the system in the horizontal direction (in lattice unites 
of T) and is given by Eq. ^ with 

P 3 (wi,u 2 ,M3) = | [P3{ui,u 2 ,u 3 ) - P^(u 1 ,u 2 ,u 3 )] 

= | [P 3 + (lii,?i2,U 3 ) +P 3 + (u 3 ,W2,Ui)] 

being the symmetric (with respect to the permutation of 
u\ and u 3 ) part of P 3 + («i, u 2 , u 3 ). The explicit expression 
for P 3 (ui, u 2 , u 3 ) which follows from Eq. I)42|) is 

P 3 °(mi,u 2 , U 3 ) = ^y'(u 3 -u 1 ) 2 (u 3 -2u 2 +u 1 ) 

+ -^[( Ul -u 2 f + (u 3 -u 2 f].{m) 

Analogous comparison of the values of P an h in two dif- 
ferent states allows to find that the free energy of a single 
type II domain wall on the background of the state (a) 
is given by 



oo 

Pdw = L x E mV(m) 

m— 1 



(47) 



whereas the interaction of two domain walls situated at 
to = toi and to = to 2 can be written as 

oo 

Pint(ni 2 - m i) = -2L X E mV(\m 2 - mi\ + m) . (48) 

m— 1 



are two lattice analogs of the derivative in the vertical 
direction suitable for a triangular lattice. It follows from 
symmetry reasons that 



<(V„u)(V±u)> = 



(50) 



The function defined by Eq. I|49() is a third order 
polynomial of variables u n>m i belonging to the stripe with 
to' = to, to + 1. Accordingly, the result of the Gaussian 
averaging of the product S° ll S% l2 will be a third order 
polynomial of the two-point correlation functions. It is 
not hard to check that the only terms which survive in 
the expression for (S mi S m2 ) are the triple products of the 
two-point correlation functions whose arguments belong 
to different stripes, whereas all other terms cancel each 
other in the result of the summation over n or as a con- 
sequence of Eq. I|5U|) . Speaking more precisely, due to the 
structure of the expression for S^, Eq. I|49l) . they are the 
triple products of the lattice analogs of the derivatives of 
such correlation functions, which for |toi — TO2I 3> 1 can 
be rather accurately calculated in the framework of the 
continuous approximation. 

When integer variables n and to are replaced by con- 
tinuous variables a; and y, both the first and the second 
term in the square brackets in Eq. 14911 should be replaced 
by the same integral 



— / (l.r ir,. 11 



2, 



where, however, u 2 x should be calculated at the values of 
y which differ from each other by one (in this subsection 
subscripts x and y designate partial derivatives with re- 
spect to the corresponding variables). This means that in 
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the framework of the continuous approximation the total 
expression for should be replaced by 

(S'.m)cont = V2K 3 J dx (ul)yUy = -\[2K 3 J dxu xx u v . 

(51) 

When writing Eq. (|51|l we have performed the integra- 
tion by parts and also the rescaling x — ► Ax (A 2 = 
1 + 2K\jK 2 — 2), which transforms the continuous ver- 
sion, 

H^l = \J Jdxdy [(2K, + K 2 )ul + K 2 u 2 y ] , 



of the harmonic Hamiltonian. 

n=m(mod2) m 

K. 



+2,m 



IT ( U »+b 



m+l 



i=±l 



to the isotropic form, 

(21 ^ 



(52) 



where 



*f eff = ^(2^1 + K 2 )K 2 = 



3/2 



Substitution of the correlation function, 

T 



G(x, y) = const 



•ln(ar + y ) l 



corresponding to the Hamiltonian i|52|l into 

1 /cont cont/ 

/OC />OC 
dxi / dx2[2G 2;:i:a ; :i: G'j (J( + 4G a , :I , y Gj / j / ] , 
-OO J —OG 

[where all derivatives of G(x, y) should be taken at x = 
x\ — x 2l y = yi — y 2 ] and subsequent integration over 
xi - x 2 give 



Vcont("l) = 7< 



T 2 



1 



where 



J \m\ — m 2 \ 7 ' 
45^ 



(53) 



= 2Kjj 15, _ 

(2ttK cS ) 3 4 1024tt 2 



0.0077 . 



Thus we have found that the quantity V(m), summa- 
tion of which over to allows to find different essential 
free energies, contains a very small numerical coefficient 
7 C and very rapidly decays with the increase of to. Ac- 
cordingly, the expressions in Eq. I|44|l and Eq. I|47|l which 
include the summation of V(m) starting from to = 1 



will be determined entirely by the first term in the sum. 
However, substitution of Eq. (|53[) into Eq. I|48|l with sub- 
sequent replacement of the summation by integration al- 
lows one to find the form of the interaction of two domain 
walls for to ^> 1, 



■Pint (m) 



— 

15 J TO 5 



F. Numerical calculations 

Note that the expression for V(m), Eq. i|53fl. which we 
have found in the framework of the continuous approxi- 
mation is positive and increases with decrease of to, that 
is when one moves out of the region of the applicability 
of the continuous approximation. Although it hardly can 
be expected that the more accurate calculation will lead 
to the change of the sign of V(m), we have checked this 
for to = 1 by going beyond the limits of the continuous 
approximation. 

The exact expression for V(m) given by Eq. I|45|) and 
Eq. (|49|l can be written as a sum over all possible pairs 
of triangles belonging to two different stripes. For to = 1 
the first term in this sum, that is the contribution which 
corresponds to the pair of adjacent triangles has the form 



_ mi 



((u„+2,m ~ Zw) 2 ) 2 ((V+w)(-V m u)) 



The averages which enter Eq. I|54|) are given by the inte- 
grals over Brillouin zone, 



((V+«)(-V-tt)) = 



dqdpW 1 {q,p)T 
_J_,{2^Y Ao(?,p) 
* dqdp W 2 {q,p)T 
„(2ir) 2 Ao( g ,p) 



where 



Wi(q,p) = 2(l-cos2<?) , 

W 2 (q,p) — (cos q — cos p) 2 — sin 2 p , 

h-o(q,p) = 2Ki(l - cos 2q) + AK 2 (l - cos g cos p) 

Numerical calculation of these integrals gives 



{{u n +2,m — Un,m) 2 ) 

((V+u)(-V-«)) 



T 

0.433— , 

-ft-2 



0.0294- 



T 
K2~ 



which after substitution into Eq. (|54|l allows one to find 
that Vb(l) = 70T 2 / J, where 70 ~ 0.0018. 

Addition to Eq. 1)54(1 of the analogous terms related to 
the pairs of triangular plaquettes which have a common 
site leads to the replacement of 70 by j 2 » 0.0032 > 0. 
The contributions from more distant pairs of plaquettes 
are much smaller and can be safely neglected. This re- 
sult confirms the positiveness of 1^(1). According to Eq. 
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(|44|l . the positiveness of V(m) for all m ensures that 



po > pa 
anh anh ' 



It follows from Eq. (|47[) that the value of the type II do- 
main wall free energy per unit length, can be then rather 
accurately estimated as 



AO) 

JDW 



72- 



J 



(55) 



In the analysis of the next section fjy^ (T) plays the role 
of the fluctuation induced effective energy of a domain 
wall. 

The comparison of the free energies of anharmonic fluc- 
tuations in the states (a) and (c) can be made following 
the same approach, but turns out to be much more cum- 
bersome for two reasons. First, in order to reduce the 
Hamiltonians of the harmonic fluctuations in these two 
states to the same form one has to apply the gauge trans- 
formation introduced in Sec. III.B. Second, there is no 
complete cancellation of the second term from Eq. (|42|l . 
which leads to the strong increase of the number of terms 
one has to take into account in the expressions for the free 
energies of fluctuations. A numerical calculation shows 
that the free energy of anharmonic fluctuations is lower 
for the state (a), which means that the free energy of the 
type I domain wall shown in Fig. 3(b) is also positive. 

The numerical constant ji characterizing this free en- 
ergy is equal to 0.0044 if one takes into account only 
the contributions from the adjacent plaquettes, whereas 
when the contributions from the pairs of plaquettes which 
have a common site are also included, one gets the value 
which is very close to 72, 

71 « 0.0033 



The temperature of the phase transition associ- 
ated with the spontaneous creation of infinite domain 
walls, T c , can be then estimated from the condition 
fuw(T c ) = 0, which can be rewritten as 



T r = 



E K 



In vT c /f£i(T ( 



(56) 



Eq. Ij56|l shows that T c is determined mainly by Ek and 
only logarithmically depends on that is on 7. 

Fig. 6(a) shows the structure of an elementary kink on 
a type I domain wall separating two different versions of 
the state (a) . This is one of the simplest neutral point- like 
excitations possible in the system. It contains only two 
vortex clusters of anomalous sizes, one with four positive 
vortices (instead of three) and another with two, so there 
is no excess vorticity associated with this defect. Note 
that any defect with only one vortex cluster of anomalous 
size will be characterized by a non-zero vorticity, so its 
energy will be logarithmically divergent. 

The distance between two neighboring positions the 
kink shown in Fig. 6(a) can occupy is equal to 2 (in lat- 
tice units of T). However, the same domain wall allows 
also for the formation of kinks of the opposite sign (ori- 




V. DISORDERING OF VORTEX PATTERN 

A. An estimate for the phase transition 
temperature 

The temperature of the phase transition associated 
with the proliferation of domain walls and the disor- 
dering of the periodic vortex pattern can be estimated 
by analyzing a more complete expression for the domain 
wall free energy, fr>w{T), which in addition to the term 
induced by anharmonicities, f^f(T) = 7T 2 /J, should 
also include the entropic term related to the formation 
of kinks, 



/dw(T) = (T) - vTexp(-E K /T) , 

where Ek oc J is the energy of a kink and v ~ 1 is the 
density (per unit length) of the positions on a domain 
wall where a kink can exist. In the case of the exactly 
solvable anisotropic Ising model j33| an analogous esti- 
mate allows one to find the transition temperature with 
the exponential accuracy. 






FIG. 6: Point-like defects with finite energies: (a) a kink on a 
type I domain wall, (b) an end point of a striped defect with 
the minimal width, (c) an intersection of a striped defect with 
a type I domain wall. 
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entation), which means that the value of v in Eq. i|56|) 
should be set equal to one. 

Numerical calculation of the kink energy has been per- 
formed by minimizing the energy of a finite lattice clus- 
ter around the kink (with the size 4L x 4L, containing 
N = 48L 2 — 10L + 1 sites inside it) with the assump- 
tion that on all sites outside of this area the values of 
the phases are exactly the same as they would be if the 
kink was infinitely far. Numerical calculation of Ek for 
L = 1, 2, 3 (L = 3 corresponds to N = 403) and extrap- 
olation of the result to L — > oo give 



J 



= 0.1037 ±0.0001 



(57) 



The simplest kink on a type II domain wall is also 
neutral, but has a more complex structure. It contains 
four vortex clusters of anomalous sizes, and, therefore, 
its energy is, roughly speaking, two times larger then the 
energy of a kink on a type I domain wall. This means 
that in the vicinity of the phase transition type I domain 
walls play a relatively more important role. Numerical 
solution of Eq. JHJ) for E K /J = 0.1037 and 7 = 0.0033 
gives 

T c /Jw 0.010. (58) 



B. Vortex pattern fluctuations in the low 
temperature phase 

At T < T c all domain walls excited as thermal fluc- 
tuations should have the form of closed loops. At tem- 
peratures well below T c the form of these loops will be 
strongly anisotropic. At T <C T c a typical defect will be 
formed by two parallel zero-energy domain walls sepa- 
rated by the minimal possible distance |2(| . Accordingly, 
the free energy of such stripe defect can be written as 



F sb (L,T) = 2Eq 



9f (0) 
Z ->DW 



(T)L 



where Eq is the energy which can be associated with each 
of the two ends of stripe defect and L is its length. Since 
in the considered problem the structure of the end point 
of a striped defect is the same as the one of the kink [see 
Fig. 6(b)], Eq is very close to Ek- 

The probability of formation of stripe defects, Psd(L), 
is, naturally, determined by their free energies, 

P SU (L) = exp[-F SD (L)/T] , 

which allows to estimate p, the fraction of the total area 
of the system covered by such defects, as 



P(T) 



T 



exp 



2Eq 
T 



==r ■ (59) 



By looking when p(T) becomes of the order of 1, a cri- 
terium is obtained which differs from Eq. i|56|l only by the 



replacement Ek — * Eq. Since Eq « Ek, this gives an ad- 
ditional support for our estimate of the phase transition 
temperature. 

However, in this approach it becomes more clear that 
the estimate we have constructed is an estimate from be- 
low. Since stripe defects are strongly anisotropic and can 
have different orientations, they have to start crossing 
each other while p is still much smaller than 1. An esti- 
mate shows that the average distance between the centers 
of stripe defects becomes comparable with their average 
length when T w |T C . Since each of such crossings costs 
an additional energy, this will decrease the rate at which 
p(T) grows with increasing temperature [as well as the 
rate at which fr>w(T) decreases]. 

The defect which is formed when a stripe defect crosses 
a type I domain wall with different orientation is shown 
in Fig. 6(c). Like the two other types of local defects 
considered above this defect is neutral and consists of 
two clusters of anomalous sizes (four and two), which 
suggests that its energy is also close to Ek- Since no 
additional energy scale is involved, one can hope that the 
effects related with such crossings will lead only to the 
appearance of some numerical factor (comparable with 
1) in the right-hand side of Eq. (|56|l . 

Note, that one also cannot exclude a possibility that 
the disordering of vortex pattern is a multistage process 
and takes place as a sequence of phase transition, the first 
of which, at T c ~ 0.01 J, is related to the appearance of 
infinite domain walls with only one orientation. 

At T — > the value of p given by Eq. exponentially 
tends to zero, which means that with the decrease of 
temperature the system becomes more and more ordered. 
Quite remarkably, this is accompanied by the divergence 
of the correlation radius of fluctuations, 



To{T) 



2/dw(T) 



(60) 



which is determined by the average length of the defect. 
In the low temperature limit r c (T) oc 1/T. 



C. Finite-size effects 

Thus, we have found that T c , the temperature of 
vortex-pattern ordering in the fully frustrated XY model 
on a dice lattice can be expected to be of the order of 
0.01 J. It has to be emphasized that at T < T c the fluc- 
tuation induced free energy of domain walls is extremely 
weak, < 10- 4 7, where 7 < 0.01 is an additional 

small parameter calculated in Sec. IV.F. 

The very low value of the ratio f^/T at T < T c 
leads to the unusual prominence of the finite-size effects 
consisting in the spontaneous formation of domain walls 
crossing the whole system. If a sample has a form of a 
stripe with a finite width, L, the probability (per unit 
length) to have a domain wall crossing the whole system 
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can be estimated as 



p(L) ~ exp 



/pw(r) 

T 



L 



exp 



L 



2r c (T) 



Vortex-pattern ordering, or, at least, any traces of such 
an ordering can be expected to be observable only when 
the average distance between such walls, r(L) = l/p(L) 
[for L < r c (T) this quantity plays the role of the effec- 
tive correlation radius induced by the finite-size effects] 
is much larger then 1 , which requires to have L>r c (T) . 

In typical systems with discrete degrees of freedom 
analyzed in statistical mechanics (the simplest example 
being the isotropic Ising model) a twofold or a three- 
fold decrease of temperature with respect to T c is usu- 
ally sufficient to obtain r c ~ 1, which allows to observe 
the ordering even in relatively small systems. However, 
in situations when a finite free energy of domain walls 
arises only from the anharmonic fluctuations, r c (T) di- 
verges not only when T — > T c , but also when T — ► 
|34|. and, therefore, the best conditions for the obser- 
vation of vortex-pattern ordering in a finite system are 
achieved at intermediate temperatures. The differentia- 
tion of fr>w(T)/T with respect to T shows that the min- 
imum of r c {T) defined by Eq. i|60[l is achieved when 



T 



\n[vE K /fZ(T)] 



(61) 



Numerical solution of Eqs. i|61[) for the same values of Ek 
and 7, and substitution of the result into Eq. H60|) show 
that the minimal value of r c is achieved when T « 0.8 T c 
and is of the order of 2 • 10 4 . Thus, the observation of 
vortex-pattern ordering requires the linear size of the sys- 
tem to be at least comparable with 10 5 . 



VI. DESTRUCTION OF PHASE COHERENCE 

Up to now we have discussed only one phase transi- 
tion related to the disordering of vortex pattern and the 
proliferation of domain walls. However, the ground states 
of uniformly frustrated XY models are characterized by 
a combination of discrete and continuous degeneracies, 
which provides possibilities for the existence of (at least) 
two different phase transitions Q . The second phase tran- 
sition is related to the vanishing of the helicity modulus, 
r(T), describing the rigidity of the system with respect 
to a phase twist. In terms of a Josephson junction array 
this phase transition corresponds to the destruction of 
superconductivity. It takes place not necessarily at the 
same temperature as the vortex-pattern disordering. 

The interaction between the discrete and continuous 
degrees of freedom in uniformly frustrated XY models 
has a non-perturbative nature and is related to the for- 
mation of fractional vortices at corners and intersections 
of domain walls 0,11 111- According to the general scheme 
proposed in Ref. 9, three scenarios are possible in a situa- 
tion when the disordering of vortex pattern takes place as 



a single phase transition (whose temperature we denote 
T c ) and not as a sequence of phase transitions (3{|. 

First, the vanishing of the helicity modulus can take 
place at T < T c , if TV, the temperature of pair dissocia- 
tion for ordinary (integer) vortices is lower than T c . The 
phase transition at T = Ty in that case has exactly the 
same nature as the Berezinskii-Kosterlitz-Thouless tran- 
sition 36] in the conventional XY model (without frus- 
tration) . Numerical simulations 0, [^} , as well as anal- 
ysis of the mutual influence between vortices and kinks 
on domain walls |13|. demonstrate that this scenario is 
realized in the fully frustrated XY models on square and 
triangular lattices. 

The vanishing of the helicity modulus can also take 
place at T > T c , but only if at T — T c the logarithmi- 
cal interaction of fractional vortices is strong enough to 
keep them bound in pairs. Note that at T < T c the con- 
finement of fractional vortices is ensured by their linear 
interaction related to a finite free energy (per unit length) 
of the domain walls which are connecting them. In such a 
case the loss of phase coherence is related to the dissocia- 
tion of pairs of logarithmically interacting fractional vor- 
tices and can be expected to take place at T = Tpv > T c , 
where T-py is the solution of the equation, 



T : 



|Q 2 r(T) 



(62) 



and Q < 1 is the topological charge of the elementary 
fractional vortex. Eq. i|62|) is nothing else but the gener- 
alization of the Nelson-Kosterlitz universal relation j^]} 
for fractional vortices |9j . This scenario is realized in the 
antifcrromagnetic XY model on a kagome lattice, where 
T r is expected to be anomalously small, T c /J ~ 1CP 4 
|34| , and also in the uniformly frustrated XY model with 
/ = 1/3 on a dice lattice, in which the vortex pattern is 
disordered at any finite temperature and becomes quasi- 
ordered only at T = pf . 

The two transitions can be expected to coincide if at 
T = T c the value of T(T) is sufficiently large to ensure 
that integer vortices are bound in pairs, but is not large 
enough to prevent from dissociation the pairs of frac- 
tional vortices. In such a case T(T) jumps to zero exactly 
at T = T C! _the ratio T/T(T) at the transition point is not 
universal [3 [(tt/2)Q 2 < T/T(T) < 7r/2L and the transi- 
tion is likely to be of the first order |g, |3!| • The results of 
numerical simulations suggest that this scenario is quite 
possibly realized on square lattice at / = 2/5 01 > as weu 
as at / = 1/8 and / = 1/10 39]. 

In the case of the fully frustrated XY model on 
a dice lattice the effective value of the helicity mod- 
ulus (properly averaged over angles) at T = is 
T = (2/3) 3 / 2 J w 0.54 J, and, therefore, at T = T c < J 
the integer vortices are strongly bound in pairs. It has 
been already mentioned in Sec. V that the kinks on both 
types of zero-energy domain walls [on the background of 
state (a)] are neutral. Therefore, the mechanism which 
forces Ty to be lower than T c in the fully frustrated mod- 
els on square and triangular lattices 1 1 .31 ] here most prob- 
ably does not work. 
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In the considered model the fractional vortices appear 
at points where two type I zero-energy domain walls cross 
each other. Fig. 7 shows an example of such a crossing. 
Note that after crossing one of the walls has to be trans- 
formed into a type II domain wall. The energy of this 
state is above the ground state energy because one of vor- 
tex clusters contains only two positive vortices instead of 
three. The accurate summation of the nominal values of 
the variables 9jk along a closed loop going around this 
cluster demonstrates the misfit of 7r/4 with respect to 
what one could expect from counting the number of pos- 
itive and negative half-vortices inside this loop. The value 
of the misfit is the same for all closed loops surrounding 
the anomalous cluster and should be compensated by a 
continuous rotation of the phase by the same angle in the 
opposite direction. This means that the cluster consist- 
ing of two positive vortices behaves itself as a fractional 
vortex whose topological charge is equal to —1/8. 

One also can construct an analogous intersection where 
one of the clusters is formed by four positive vortices in- 
stead of three. The topological charge of such a defect 
will be equal to +1/8. It is clear that when the cluster 
of an anomalous size consists of negative vortices (in- 
stead of positive), the sign of the topological charge is 
reversed. The topological charges of more complex inter- 
sections (containing, for example, the clusters of five or 
more vortices or several clusters of anomalous sizes) will 
all be the multiples of 1/8. Note that the excess vorticity 
which can be associated with a vortex cluster depends 
not only on its size but also on the shape. For exam- 
ple, when a cluster consists of three vortices, but has the 
shape of a hexagon, the topological charge which has to 
be associated with it is equal to ±3/8. 

The only possibility to have domain walls crossings in 
an equilibrium infinite system well below the temperature 
of vortex-pattern ordering is related to crossing of stripe 
defects discussed in Sec. V.B. Fig. 6(c) shows an exam- 
ple of the intersection of a stripe defect with a domain 
wall where the two fractional vortices have the opposite 
topological charges, which makes the energy of such an 
object finite. Although it is possible also to construct an 
example in which the topological charges of the two in- 
tersections will be the same, the total topological charge 




FIG. 7: An example of a fractional vortex 



of any finite size defect (for example, formed by several 
intersecting stripe defects) has to be integer pi Hi - There- 
fore, at low temperatures the fractional vortices can be 
present only in the form of bound pairs, whose size is 
restricted more by the available separations between do- 
main walls in stripe defects rather then by the logarith- 
mic interaction of these objects (which is 64 times weaker 
than the interaction of ordinary vortices). 

With increase of temperature the average separation 
between the domain walls forming a stripe defect in- 
creases, which allows larger separations between the frac- 
tional vortices of opposite sign. Above the temperature 
of vortex-pattern disordering, the restrictions for the dis- 
tances between fractional vortices related with the sep- 
arations between domain walls in stripe defects will no 
longer exist. It looks rather likely |9j that in such a situa- 
tion one can consider fractional vortices (at the scales 
that are large in comparison with the correlation ra- 
dius) as really logarithmically interacting objects. The 
same approach may be also applicable to a finite system 
at T < T c if its linear size does not exceed the size- 
dependent correlation radius r c {L) introduced in Sec. 
V.C. In both cases one can speculate about the possi- 
bility of a phase transition related with the unbinding of 
neutral pairs formed by logarithmically interacting vor- 
tex clusters of anomalous sizes. 

Substitution of Q — 1/8 into Eq. I|62l) allows one to 
find that the temperature of this phase transition can be 
estimated as 

Tfv » ^r » 0.013 J . 

This is an estimate from above which neglects the renor- 
malization of T by thermal fluctuations. Comparison with 
the estimate T c > 0.01 J obtained in Sec. VI suggests 
that in the fully frustrated XY model on an infinite dice 
lattice the destruction of phase coherence will be trig- 
gered by the disordering of vortex pattern, which can be 
expected to occur at a temperature where the logarith- 
mic interaction of fractional vortices at T c is too weak to 
keep them bound in pairs. 

On the other hand, in a situation when the size of 
the system is insufficient to exclude the specific finite- 
size effects leading to the disordering of vortex pattern 
at any temperature (see Sec. V.C), one still can discuss 
the possibility of a phase transition (slightly smeared by 
the finite-size effects), in which the loss of phase coher- 
ence will occur as the result of the unbinding of fractional 
vortices at T = Tpv ~ 0.01 J. In numerical simulations 
this phase transition can be observed by analyzing if vor- 
tex clusters of anomalous sizes are bound in neutral pairs 
or not. However, at T ~ 0.01 J one should be specially 
attentive about checking if the time of simulation is suffi- 
cient for the equilibration of the vortex subsystem, which 
at such temperatures will require much longer times than 
the equilibration of the spin-wave subsystem. It is not 
clear if in numerical simulations of Ref. [lj (discussed in 
a more detail in Sec. VIII) this condition was really sat- 
isfied. 
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VII. MAGNETIC EFFECTS 

It has been already mentioned in the Introduction that 
the main interest to the uniformly frustrated XY mod- 
els has appeared in relation to their application for the 
description of Josephson junction arrays in external mag- 
netic field. Since we have found that in the case of the 
fully frustrated model on a dice lattice the order-from- 
disorder mechanism for the removal of an accidental de- 
generacy is extremely inefficient, in a physical situation 
one should also take into account other possible mecha- 
nisms. In the case of a proximity coupled array, the most 
important of them is related to the magnetic interactions 
of currents IM EH . 

When the proper magnetic fields of currents in the ar- 
ray are taken into account, the Hamiltonian of the frus- 
trated XY model should be replaced 0>E3 by 

H = - COs(6>jfc - Ojfe) + -Emagn (63) 

W 

where 9jk = tfk — — Ajk includes only the contribu- 
tion from the external magnetic field, defined by Eq. J5J, 
whereas the analogous contribution from the currents is 
denoted ajk- The second term in Eq. l|6lfl) . 



E 



magn — „ ^ afi 

a,/3 



(64) 



is the energy of the current induced (screening) magnetic 
fields expressed in terms of 



(65) 



the magnetic fluxes of these fields through the plaque- 
ttes of the array (denoted by Greek letters), L a p being 
the matrix of mutual inductances |44j between the pla- 
qucttes. 

The values of the variables ajk should be found from 
the minimization of H. The result of the variation of Eq. 
()63|l with respect to ajk can be rewritten as 
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(66) 



where Ip is so-called mesh current |44j, |45j which can 
be associated with the plaquette /3. The currents in the 
junctions, 

Ijk = h sin(8jk - ajk) , 

are given by the difference of mesh currents in the two 
plaquettes, ajk and a'j k , sharing the bond (jk), 



Ijk — Ia jk I a 



(67) 



The substitution of Eq. ifHHf into Eq. (|6*4*j) allows to 
rewrite it in the more familiar form as 



E 



1 



\ " 



magn 



2 ^ 



In the regime of weak screening one can assume that 
the values of the currents are not affected by their 
magnetic fields and calculate -Emagn replacing Ijf. by 
ijJP = Iosm8jk- Nonetheless, the calculation of the first 
term in Eq. I|63|l with the same accuracy requires to take 
into account its dependence of ajk- The contribution to 
this term which is linear in ajk has a form 

(0), 



m 

and with the help of Eqs. I|64|) - (|67|l can be shown to 
be equal to — 2£ , magn , where -E mag n is calculated for the 

"bare" values of currents, 1^, . Therefore, in the regime 
of weak screening, H niagn , the total magnetic correction 
to the Hamiltonian of a Josephson junction array is equal 
to -Emagn, but has the opposite sign, 



-E, 



magn • 



The comparison of -E magn in different periodic states in 
a fully frustrated superconducting wire network with the 
dice lattice geometry has been recently made in Ref. l24l 
Since in terms of the gauge-invariant phase differences 8jk 
the structure of these states is in one-to-one correspon- 
dence with the ground states of the fully frustrated XY 
model with the same geometry, the same results are also 
applicable to an array formed by weakly coupled super- 
conducting islands. It follows from the results of Ref. 00] 
that in arrays the expression for -Emagn (normalized per a 
sixfold coordinated site) in the regime of weak screening 
can be written as 

e J 2 

Bmara = — - (68) 



where 



3E$ 



* 2 o 
in 2 a 



is the characteristic energy scale which depends on a, the 
lattice constant of the array. The numerical coefficient e 
in Eq. I)68|l is determined by the structure of vortex pat- 
tern in the considered ground state and can be expressed 
as a linear combination of the coefficients Xi > defining 
the mutual inductances, Li = — A^a, between different 
pairs of plaquettes on a dice lattice. Here i = 1 corre- 
sponds to the nearest neighbors, i = 2 to the next-to- 
nearest neighbors, etc.. 

The main conclusion of Ref. |3 is thsit, in the fully 
frustrated system with the dice lattice geometry the co- 
efficient e is the largest in the state (c). For example, the 
numerical coefficient in the expression, 

I 2 

err _ rra rrc _ " 

magn magn magn ' 

for the difference between the magnetic energies of the 
states (a) and (c), 

e c - e a 1 
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is positive. 

For a = 8 (29| ps 0.98 • 10 4 K, which shows 
that at T ~ J the differences between the magnetic en- 
ergies are extremely small, SH masn /T < 10~ 4 , and are 
even smaller than the differences between the free ener- 
gies of anharmonic fluctuations. In this estimate we have 
assumed T < 10 K. 

It should be emphasized that in proximity coupled ar- 
rays arrays the coupling constant J has a strong tempera- 
ture dependence in a wide interval of temperatures, so the 
decrease of the dimensionless temperature t — T/J(T) 
with the decrease of real temperature T is much more 
strongly influenced by the growth of J(T) than by the 
decrease of T. Roughly speaking, in the experiments of 
Ref. |29J the decrease of T by 1 K corresponds to the de- 
crease of r by one order of magnitude. 

This suggests that the importance of magnetic effects 
rapidly grows with the decrease of r. Comparison of 
SH magn with SF anh = 7T 2 / J shows that the magnetic 
energy becomes dominant when 

T < T mag n = — — P3 0.30 , 

V7 

where we have put T = 5K. Thus at r < 0.1, where 
the vortex-pattern ordering can be expected to occur, 
one can take into account only the magnetic energies of 
different states (or domain walls), completely neglecting 
the free energy of anharmonic fluctuations. Therefore, 
the structure of the periodic vortex pattern in the low 
temperature phase of a proximity coupled array should 
be of the type (c). 

The energy of domain walls separating different ver- 
sions of state (c) from each other, Ejjw, will be close to 
SH ma ,g n , which will make the vortex pattern more robust 
with respect to thermal fluctuations in comparison with 
the XY model. However, any quantitative conclusions 
about the temperature of the phase transition(s) related 
to vortex-pattern disordering in this situation are im- 
possible without the detailed analysis of the energies and 
other properties of the topological excitations in state (c) 
[analogous to the one performed in Sec. V for state (a)], 
which goes beyond the limits of this article. Nonethe- 
less, some increase of the region of stability of the low- 
temperature phase with the ordered vortex pattern is in- 
evitable, which allows to conclude that the scenario in 
which the unbinding of fractional vortices with topolog- 
ical charges ±1/8 takes place as an independent phase 
transition is impossible. 

Inclusion of the magnetic interactions into analysis im- 
proves also the situation with respect to the finite-size 
effects. For £"dw ~ J 2 /-E$ the probability to have a do- 
main wall crossing the whole system (of the width L) 
becomes negligible for 

L > L c ~ — r 2 , 

which at r <C 1 is a much more mild restriction than the 
one obtained in Sec. V for the XY model. 



On the other hand, it is known that the magnetic inter- 
actions of currents, or, in other words, screening effects, 
lead to the increase of barriers a vortex has to overcome 
when moving between the plaquettes of an array |45| . 
Thus, although the growth of screening effects with a de- 
crease in temperature improves the stability of a vortex- 
pattern ordering, it simultaneously leads to the further 
increase of (already exponentially large) times required 
for the relaxation of the system to equilibrium. This may 
be one of the reasons why the finite-frequency experi- 
ments of Ref. have not demonstrated any signs of a 
phase transition at / = 1/2. 

It should also be noted that in the case when the loss 
of phase coherence is related to dissociation of fractional- 
vortex pairs, the universal relation [46] for the value of 
the two-dimensional magnetic penetration depth |47l ]. A, 
at the transition temperature can be rewritten (in our 
notation) as 

A(Jfv) = Q 2 E* 
a 8 Tpv 

As a consequence, the smearing of the phase transition 
due to the finiteness of A in the case when it is related 
to unbinding of fractional vortices should be more pro- 
nounced then in the case of integer vortices. 

VIII. CONCLUSION 

In the present work we have investigated the fully frus- 
trated XY model on a dice lattice and have demonstrated 
that the energy of this system is minimized in the highly 
degenerate family of states (described in Sec. II. C), in 
which the vortices of the same sign are grouped into clus- 
ters of three. The accidental degeneracy of these states 
can be described in terms of the formation of a network of 
intersecting zero-energy domain walls, whereas the resid- 
ual entropy related to this degeneracy is proportional to 
the linear size of the system, as in the case of a honey- 
comb lattice 

The central result of this article consists in determin- 
ing the structure of the periodic vortex pattern which is 
selected at low temperatures by thermal fluctuations. It- 
is shown in Fig. 3(a). However, this effect is rather weak, 
being induced only by the anharmonic fluctuations. As a 
consequence of a hidden gauge symmetry, the free energy 
of the harmonic fluctuations turns out to be the same for 
all ground states. The same conclusion is applicable also 
to quantum generalizations of the model both at a fi- 
nite and at zero temperature, when one should speak of 
zero-point fluctuations. 

The destruction of the periodic vortex pattern with the 
increase of temperature is related to the proliferation of 
domain walls. The dimensionless temperature, r = T/J, 
at which the corresponding phase transition can be ex- 
pected to take place can be estimated as r c ~ 0.01, where 
one factor of 0.1 is related to the smallness of the energy 
of particular point-like defects on domain walls, and the 
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other (which is of the logarithmic origin) to the extreme 
smallness of the fluctuation induced free energy of zero- 
energy domains walls. The analysis of possible scenarios 
suggests that the loss of phase coherence in the consid- 
ered system can be expected to take place at the same 
temperature as the disordering of vortex pattern. 

The extreme smallness of the fluctuation induced free 
energy of domain walls will manifest itself also in the 
huge prominence of the finite-size effects consisting in 
the appearance of domain walls crossing the whole sys- 
tem and leading to the disordering of vortex pattern even 
at r < t c . As a consequence, even at the "optimal" tem- 
perature, t w 0.8t c , the linear size of the system required 
to observe the vortex-pattern ordering should be much 
larger than r™ m « 2 • 10 4 . However, in smaller systems 
one can still discuss a possibility for the observation of 
a phase transition related to the loss of phase coherence 
and consisting in unbinding of pairs of logarithmically 
interacting vortex clusters of anomalous sizes, which can 
be expected to happen at tfv ~ 0.01. 

Our conclusions are consistent with the results of the 
numerical simulations of the same model by Cataudella 
and Fazio |14j . who have investigated the temperatures 
down to t = 0.01 and have found no signs of vortex- 
pattern ordering. The results presented in Sec. V demon- 
strate that even if the lowest temperatures investigated 
in Ref. were indeed below r c , the size of the system 
was definitely not sufficient for the observation of such 
an ordering. In the notation based on counting only the 
sixfold coordinated sites (used in this work), the largest 
size of the system analyzed in Ref. ^3 corresponds to 
56 x 42, whereas the majority of the data has been taken 
at 24 x 18. It seems rather likely that analogous simula- 
tions of a 10 5 x 10 5 system, which may be required for 
the observation of vortex-pattern ordering in the fully 
frustrated XY model on a dice lattice, can become pos- 
sible only with a further development of computational 
abilities. 

Another result of the numerical simulations of Ref. 
fill is related to the dynamical properties of the system 
(which were not discussed in this article) and consists in 
finding below r* s» 0.06 the signs of a glassy behavior. 
Namely, the relaxation of the total energy at t < r» be- 
comes logarithmic in time in contrast to the exponential 
relaxation at t > r,, whereas the behavior of the auto- 
correlation function of the variables 6jk demonstrates a 
dependence on the waiting time. 

One can certainly make an attempt to relate this obser- 
vation to the specific features of the considered system. 
For the temperatures and system sizes analyzed in Ref.ll4l 
one can safely neglect the fluctuation induced free energy 
of zero-energy domain walls and treat them as objects 
with zero free energy. It seems rather likely then that the 
typical state obtained after cooling down the system from 
t ~ 1 to t <C Ek I J can be described as a network formed 
by zero-energy domain walls, which contains a large num- 
ber of point-like defects, such as kinks and intersections, 
which cost an additional energy. We know that the inter- 



sections with zero energy are also possible [see Fig. 3(f)], 
but one can expect that the majority of the intersections 
formed during cooling down from a disordered state will 
not have the optimal structure necessary for that. 

The glassy behavior can be then expected from the 
necessity of the disentanglement of this domain wall net- 
work with the decrease of temperature. For different 
temperatures above r c , the equilibrium concentration of 
point-like defects should be different. However, in con- 
trast to kinks on a domain wall whose number can be 
changed due to annihilation of two kinks of opposite signs 
(which may be a relatively fast process), to change the 
number of the intersections of domain walls one has to 
change the number of these walls, which will require much 
longer times than the annihilation of kinks. In particu- 
lar, it seems rather likely that the processes related to 
the annihilation of zero-energy domain walls crossing the 
whole system will be characterized by relaxation times 
diverging with the system size, providing thus a source 
for a genuine glass-like behavior. 

On the other hand, even if the times characterizing 
the relaxation of vortex pattern are not divergent in the 
thermodynamic limit and are restricted only by the bar- 
riers whose typical height remains of the order of J, this 
already can be the source for a glassy-like behavior (as- 
sociated with a wide distribution of such barriers) at 
t < 1 in a wide interval interval of times. For example, 
exp(l/0.05) ~ 10 9 , whereas in Ref. Q the anomalous re- 
laxation has been studied only at much shorter times. In 
such a case one can argue that the observation of a glassy- 
like behavior in the considered XY model is possible due 
to a combination of three factors, namely, (i) the exis- 
tence of zero-energy domain walls, (ii) the special inef- 
fectiveness of the order-from-disorder mechanism for the 
removal of an accidental degeneracy and (iii) the anoma- 
lously low transition temperature (t c ~ 0.01). However, 
the choice between the two scenarios (genuine glass vs. 
the dynamical freezing of vortex relaxation) should be 
made on the basis of a much more detailed analysis of 
the domain- walls disentanglement. 

In the experimental situation, the magnetic interac- 
tions of currents in a Josephson junction array will be 
of greater importance for the stabilization of a particular 
vortex pattern then the anharmonic fluctuations. This 
mechanism leads to the selection of the pattern shown 
in Fig. 3(c) and makes the periodic vortex pattern less 
vulnerable with respect to fluctuations. 
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